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ABSTRACT 

We show that supersonic MHD turbulence yields a star formation rate (SFR) as low as observed in 
molecular clouds (MCs), for characteristic values of the free-fall time divided by the dynamical time, 
tg/tdyn, the alfvenic Mach number, Ai a , and the sonic Mach number, M. s . Using a very large set 
of deep adaptive-mesh-refinement simulations, we quantify the dependence of the SFR per free-fall 
time, eg, on the above parameters. Our main results are: i) eg decreases exponentially with increasing 
iff/idyn, but is insensitive to changes in Ai s , for constant values of iff/^dyn and M.&. ii) Decreasing 
values of VW a (stronger magnetic fields) reduce eg, but only to a point, beyond which eg increases 
with a further decrease of iii) For values of M. & characteristic of star-forming regions, eg varies 
with A4 a by less than a factor of two. We propose a simple star-formation law, based on the empirical 
fit to the minimum eg, and depending only on tg/td yn '- e ff ~ e wind exp(— 1.6 tg /tdyn)- Because it 
only depends on the mean gas density and rms velocity, this law is straightforward to implement in 
simulations and analytical models of galaxy formation and evolution. 

Subject headings: ISM: kinematics and dynamics — stars: formation — magnetohydrodynamics — 
turbulence 



1. INTRODUCTION 

Stars are the result of the gravitational collapse of cold 
gas, but on large scales the conversion of cold gas into 
stars takes much longer than a collapse time. Early 
explanations for this low star formation rate (SFR) in- 
volved the magnetic support against the collapse of gi- 
ant molecular clouds (GMCs), while the more recent sce- 
nario of turbulent fragmentation r elies on the support of 
GMC s by supersonic turbulence (Mac Low & Klessen 
2004 1 . However, numerical simulations of star formation 



m turbulent clouds have generally failed to yield SFRs 
as low as observed valueaM (see § 4). 

This work shows that supersonic MHD turbulence re- 
sults in a SFR as low as observed in GMCs, for charac- 
teristic values of the free-fall time divided by the dynam- 
ical time, tg/tdyn, the alfvenic Mach number, Ai a , and 
the sonic Mach number, A4 S (defined in § 2). Using a 
very large set of adaptive-mesh-refinement (AMR) simu- 
lations of driven MHD turbulence, including self-gravity 
and sink particles, we quantify the dependence of the 
SFR per free-fall time, eg, on tg/t^yn, confirming the re- 
sults of our previous uniform grid simulations in the case 
of a very weak mean magnetic field. 

The SFR per free-fall time is the efficiency factor of a 



1 After submitting this Letter, we became aware of a closely 
related paper by Fcdcrrath and Klessen, submitted to the Astro- 
physical Journal (arXiv:1209.2856), presenting a large number of 
star formation simulations that are complementary to ours. 



theoretical Schmidt-Kennicut law, 

Pstars = CffPgasAff, (1) 

often used in models and simulations of star formation 
in galaxies. We find that eg decreases exponentially with 
increasing tg/td yn , has no dependence on M. s alone, and 
varies by less than a factor of two with variations in M. a 
within reasonable values for star-forming regions. We 
therefore propose a new empirical law of star formation 
that depends only on tg/t^ yn . 

2. NUMERICAL SIMULATIONS 

The si mulations are carried out with the AMR code 
Ramses ( Teyssier|2002 ), modified to include random tur- 
bulence driving and sink particles, and optimized for 
supercomputers with large multi-cor e nodes through an 
OpenMP/MPI hyb rid layout. As in |Padoan fc Nordlund" 
( |2002| |2004| |2011[ ), we adopt periodic boundary condi- 
tions, isothermal equation of state, and solenoidal ran- 
dom forcing in Fourier space at wavenumbers 1 < fc < 2 
(fc = 1 corresponds to the computational box size). We 
choose a solenoidal force to make sure that dense, star- 
forming structures are naturally generated in the turbu- 
lent flow, rather than directly imposed by the external 
force. 

The simulations are first run without self-gravity for 
approximately five dynamical times, starting with uni- 
form density and magnetic fields, and random veloc- 
ity field with power only at wavenumbers 1 < fc < 2. 
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TABLE 1 

Non-dimensional parameters and eg of the simulations. 
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The driving force keeps the rms sonic Mach number, 
A^s = ov,3d/cs, at the approximate values of either 10 
or 20, characteristic of GMCs on scales of order 10 pc 
(cv,3D is the three-dimensional rms velocity, and c s is 
the isothermal speed of sound). For the simulations with 
A4 S w 10, we have chosen two different values of the mean 
magnetic field, corresponding to rms alfvenic Mach num- 
bers of A4 a ~ 33 and 5, where A^ a = &v,3d/v & , and w a is 
the Alfven velocity corresponding to the mean magnetic 
field and the mean density. For the larger Mach number, 
M s ~ 20, we have adopted three different values of the 
mean magnetic field, corresponding to M a « 20, 5, and 
1.25. 



In the supersonic regime, th e turbulent dynamo is sup- 
pressed (Haugen et al. 2004), and the magnetic energy 
growth saturates at approximately 2% of the kinetic en- 
ergy of the turbulence, when the mean magnetic field is 



negligible and d oes not characterize the system (Feder 
[rath et al. 2011 ). On the other hand, if the mean mag- 



netic field is not negligible, as in our simulations and in 
real GMCs, it is readily amplified, well beyond the dy- 
namo saturation level, by stretching and compression in 
the turbulent flow, and thus its saturated rms value de- 
pe nds directly on its m ean value, as found in simu lations 
of|Kritsuk et al.|(|2009 ), and in equation (20) of Padoan 
&: IN ordlund 12011 ). 'therefore, the value of Ma,, defined 
with the mean magnetic field, characterizes the saturated 
regime and can be taken as a fundamental parameter of 
the system. 

The final snapshots of these five turbulence simulations 
are used as initial conditions for the star-formation sim- 
ulations with self-gravity and sink particles. Each initial 
condition is used for 5 to 7 simulations with different nu- 
merical values of the gravitational constant, correspond- 
ing to different values of tff/^dyn, resulting in 28 star- 
formation simulations (plus 6 larger runs used to test 
numerical convergence, and 11 smaller runs to test both 
convergence and variations of eg with initial conditions 
- a total of 45 star- formation simulations). The dynam- 
ical time is defined as idyn = £/2ov.3D, where L is the 
size of the computational domain, and tg = -\/37r/32G/9n, 
where po is the mean density in the computational vol- 
ume. The simulations were run until the star formation 
efficiency value of SFE = M±/M w 0.2-0.3, where M* is 
the mass in the sink particles, and M the total mass (gas 
plus sinks), except the simulations with very low SFR, 
that were run until SFEsa 0.05-0.1. 

The five initial turbulence simulations have a root grid 
of 128 3 computational cells and 5 AMR levels, thus a 
maximum resolution equivalent to 4,096 3 cells. In all 
runs, the first AMR level is created if the density is 2.5 
times larger than the mean. Subsequent AMR levels are 
created when the gas density increases by a factor of 4 
rela tive to the previous level. I n con trast to |Kritsuk et] 



al. (2006) and Schmidt et al. (2009), we do not refine 



based on velocity or pressure gradients. However, we 
verified that the tails of the density PDFs are sufficiently 
converged in the turbulence simulations used as initial 
conditions (details will be given elsewhere). When self- 
gravity is included, the root grid is reduced to 64 3 cells, 
while the number of AMR levels is increased from 5 to 8, 
to achieve a maximum resolution equivalent to 16,384 3 
cells. 

Sink particles are created in cells where the gas den- 
sity is larger than 10 times the mean density (much 
larger than the maximum density reached in the initial- 
ization runs without self-gravity), if the following addi- 
tional conditions are met at the cell location: i) The 
gravitational potential has a local minimum value, ii) the 
three-dimensional velocity divergence is negative, and iii) 
no other previously created sink particle is present within 
an exclusion radius, r cxc i (r cxc i = 17.5Ax in these sim- 
ulations, where Ax is the highest-reso lution cell size). 
These conditions are similar to those in Fcderrat h et al.l 
( |2010| . We have verified that they avoid the creation 
ot spurious sink particles in regions where the gas is not 
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Fig. 1. — Projected density, N co \, and sink particles from the 
run 128M10B2G50, at a time when SFErs 0.05. The color scale 
is linear in TV 0,2 , and the tick values in the color bar give N co \ in 
numerical units, with mean gas density and computational box size 
equal to unity. 

collapsing. 

The numerical parameters characterizing the simula- 
tions are given in Table 1. The simulation names con- 
tain the rms sonic Mach number, and the numerical val- 
ues of the mean magnetic field and of the gravitational 
constant. In Table 1, we give three non-dimensional pa- 
rameters, A4 S , .M a , and tg/td yn , that fully characterize 
the simulations. The ratio fff/tdyn expresses the rela- 
tive strength of gravity and turbulence, like the virial 
parameter, a V h-, used in previous works. Although for- 
mally equivalent to a v i r , for example through the usual 
expressions valid for a uniform sphere, we prefer to use 
ifr/idyn as the third parameter, in line with the choice of 
the Mach numbers (ratios of crossing times as well) as 
the first two parameters. 

In Table 1, we also give the ratio between the mini- 
mum Jeans length at each level, Lj^ (corresponding to 
the highest density before further refinement), and the 
computational cell size, Axi, at that level. Due to the 
refinement by a factor of two in size for every increase by 
a factor of four in density, the minimum number of cells 
per Jeans length, ijj/Aa;i, is the same for every level i 
in a given simulation. The last column of Table 1 shows 
the values of eg obtained from the simulations. 

At the bottom of Table 1, we also list the simulations 
used to test numerical convergence and to study the vari- 
ations of eg with initial conditions. The names of these 
extra runs begin with their root grid size, either 128 3 or 
32 3 computational cells. Projected density and sink par- 
ticles from one of these runs, 128M10B2G50, are shown 
in Figure 1, at a time when SFEw 0.05. 

3. RESULTS 

Figure 2 shows the time dependence of the SFE for the 
two sets of simulations with Ai s ~ 20 and A4 a ~ 5 and 
1.25. In some runs with high eg, the SFR tends to grow 
with time. To extract a representative value, we measure 
eg as the slope of the least-squares fit of the SFE versus 
t/tg, in the interval 0.03 < SFE < 0.2, as in |Padoan fc 
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Fig. 2. — SFE versus time for the two sets of simulations with 
M s w 20 and Ma, « 5 and 1.25. The values of t B /t dyn are 0.77, 
1.09, 1.54, 2.18, and 3.09, from top (high SFE) to bottom. 



Nordlund (2011 ). The simulations with the lowest values 
of eg were not integrated up to SFE=0.2. However, they 
do not exhibit a systematic growth of the SFR with time, 
and they are run for many free-fall times (up to 9 in some 
cases), so their average eg is nonetheless well defined. 

The values of eg from all simulations are plotted in Fig- 
ure 3, as a function of tg/td yn - A striking result, stressed 
by the linear-log scale of Figure 3, is the approximately 
exponential dependence of eg on tg /td yn (for fixed values 
of A4 S and A4 a ), showing that supersonic MHD turbu- 
lence is very effective at slowing down the star formation 
process, which seems to be controlled primarily by the ra- 
tio of the characteristic times of gravity and turbulence. 

Another important result is the lack of a dependence 
of eg on Ai s , for fixed values of Aia, and tg/td yn - Fur- 
thermore, the dependence on A4 a is not very strong, con- 
sidering the very large range of A4 a values of the simu- 
lations. This is in part due to the complex dependence 
of eg onM a . Down to M a ~ 5, eg decreases with de- 
creasing VW a (increasing magnetic field strength), while 
for further decreasing Ai a , eg starts to increase. 

Figure 3 shows that a simple exponential law, where 
eg is only controlled by the ratio tg/tdyn > could describe 
all simulation results within a factor of two. For even 
weaker magnetic fields than considered here, eg would 
certainly be larger, as it should approac h the values found 
in simulati ons without magnetic fields (Padoan & Nord 
lund 2011). However, observations in star torming re 
gions are consistent with M a sa 1 — 10 j Padoan fc Nord 



lund 



[19991 ILunttila et ah] [20081 [2009l THey. 

) , so it is likely that eg approaches it 



er & Brunt 



20121, so it is likely that eg approaches its minimum 
value, which is approximately fit by a simple exponential 
law, eg R3 exp(— 1.6ifi/idyn), shown as a dashed line in 
Figure 3. Wc therefore propose a new empirical law of 
star formation, based on this fit to the minimum eg: 

£g ~ ^wind exp(-1.6t ff /£ d yn), (2) 

where e W i n d ~ 0.5 accounts for mass loss through jets 
and winds, during the formation of a star. This law only 
depends on the mean gas density and the rms velocity 
of a star-forming region, so it is easily implemented in 
analytical models and simulations of galaxy formation 
and evolution. 
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Fig. 3. — SFR per free-fall time, eg, versus tff/tdyn- The symbols 
for each series of runs, where only the strength of gravity is changed 
and Ai B and Aia are kept constant, are connected by a line, to 
better distinguish each series. The two error bars give the mean 
of eg, plus and minus the rms values, for each group of five 32 3 - 
root-grid runs with identical parameters (M B ~ 10 and Ai a ~ 5), 
but different initial conditions. The dashed line is an approximate 
exponential fit to the minimum value of ejf versus iff /t^ yn . 

Equation pi) is valid as long as the fraction of the sys- 
tem mass, Mo, available to form stars is larger than the 
characteristic stellar mass, m*, which we can express 
as effMo/m* > 1. If eg predicted by equation (f2j) is 
too small to satisfy this condition, then eg = 0, simply 
because the total system mass is too small. Adopting 
the characteristic stell ar mass predicted by our turbulent 
fragmentat ion model ( Padoan fc Nordlund|2002 |Padoan 
et al. 2007 ) , 3Mbe,o/A1s, which corresponds to ap- 



proximately three times the mean Bonnor-Ebert mass, 
Mbe,0) evaluated for an external pressure equal to the 
ram pressure of the turbulence, and using the minimum 
eg derived here, we find the condition can be expressed 
as egM /m* « eg(tg/t dyn )~ 3 Aff > 1. 

As shown by Table 1, the weakest gravity we have sim- 
ulated corresponds to tg/tdyn rj 3.09, giving eg-Mo/m* rj 
2.4 and 38.6 for the runs with Ai s ~ 10 and 20 respec- 
tively. In good agreement with these values, in each of 
the A^ s ~ 10 runs with tg/tdyn ~ 3.09, we get only a sin- 
gle star, so the SFR is practically zero, while in the three 
M. a ~ 20 runs with tff/tdyn ~ 3.09 we get approximately 
10-20 stars per free-fall time. 

To test the numerical convergence of eg, we have re- 
run three of the simulations with Ai s ps 10 and A4 a f=s 5, 
using both twice smaller and twice larger root grid size 
(32 3 and 128 3 computational cells), and the same num- 
ber of AMR levels, corresponding to maximum resolu- 
tions equivalent to 8,192 3 cells and 32,768 3 cells. These 
three larger runs were then rerun again using a 10 times 
larger value of the threshold density above which sink 
particles are created (10 6 times the mean density, in- 
stead of 10 5 ). As shown in Table 1, the values of eg in 
the runs with tff/tdyn = 0.39 and 1.09 are practically 
independent of numerical resolution, showing that the 
corre ct SFR can be derive d even when Truelove's condi- 
tion ( Truelove et al.|1997 ) is violated, and that the value 
of 10° times the mean density (adopted in the main pa- 
rameter study), above which sink particles are allowed 
to form, is sufficiently high to avoid dense regions that 



are not collapsing. 

In the case of tg/tdyn = 1-54, eg is approximately 20% 
lower in the 128 3 -root-grid simulations than in the cor- 
responding 64 3 -root-grid run. However, this low SFR is 
obtained from an average over many dynamical times. 
During the first dynamical time, the SFR is virtually in- 
dependent of resolution. The corresponding simulation 
with 32 3 root-grid size has the same value of eg as the 
largest simulations, showing that the difference from the 
64 3 -root-grid run is not due to lack of convergence, but 
more likely to random differences between runs after a 
long time, as expected in chaotic systems. The values 
of eg from the 128 3 -root-grid simulations can be seen in 
Figure 3. 

To investigate the effect of the initial conditions, we 
have carried out 10 simulations with twice smaller root 
grid size (32 3 computational cells), and the same num- 
ber of AMR levels, five of them with tff/td yn = 1.09, 
and the other five with tg/tdyn = 1-54, starting from 
different initial conditions of fully-developed turbulence, 
with M s ~ 10 and M a ~ 5. As shown in Table 1, the 
maximum difference in eg from such small samples is ap- 
proximately a factor of two. The mean values of eg, plus 
and minus the rms values, are plotted as error bars in 
Figure 3. 

4. DISCUSSION 

The values of eg predicted by equation ^ ar e consis- 
tent with those fou nd in GMCs. As discussed in IPadoanl 
fc Nordlund| ( |2011[ ), accounting for the full uncertainty 
in the life time of the Class II phase (2 ± 1 Myr), the 
results of |Evans et~al] ([20 09)) implies eg = .02-0. 12. 
Using the CMC sample by fHeyer et al.| ( |2009| ), we de- 
rive tg/tdyn ~ 1.4 ± 0.5 (l-o), resulting m eg = 0.02- 
0.12 based on our equation pj), with e w j n d ~ 0.5, in nice 



agreement with the values derived by |Evans et al. (20091 
(or tff/tdyn ~ 1.0 ± 0.4 and eg = 0.06-0.19 it we double 
the cloud masses). Similar values are found also for the 
Orion Nebula Cluster, e g = 0.03-0.09, the obje ct with 
the smallest error bars in |Krumholz fc Tan fl2007[). Based 
on the values of tg/tdyn m the sample ot Murray (2011 1 
we would predict an average value of eg = 0.13 ± 0.0'/, 
while, based on the free-free emission, Murray finds an 
average value of eg = 0.14 ± 0.14. 

The AMR simulations with M s « 10 and M a 33 
presented here a re consistent with our pre vious uniform 
grid simulations (Padoan & Nordlund 2011 ) with M s ~ 9 
and M a ~ 30 (the values ot a vir given m the plot of 
Figure 2 are derived from tg/tdyn based on the case of 
a uniform sphere, while t hose reported in Figure 7 of 



Padoan & Nordlund (2011 ) were based on the total mass 
in the computational box, and should be multiplied by a 
factor of 2 for the c omparison) . 



The SFR model of |Padoan fc Nordlund| j20lT] ) adopted 
the simplifying assumption that, in the MUD case, the 
effective postshock ratio of gas to magnetic pressure, /3, 
and the density pdf do not depend on the mean mag- 
netic field strength, for reasonable values of A^a- The 
complex, but weak dependence of eg on Af a found in 
the simulations justifies that approximation. However, 
by neglecting the dependence of f3 on both M a and M s , 
the model yields a direct dependence on Ai s , not found in 
our AMR simulations. Equation (20) of Padoan & Nord 



lund ( |2011[ ), relating (3 to M s and /3q, the ratio ot gas 



■5 



to magnetic pressure computed with the average mag- 
netic field strength, oc 0n M7 1 , was derived under 



the assumption of flux freezing, neglecting the postshock 
thermal pressure a nd dynamical alignment of m agnetic 
and velocity fields. Padoan & Nordlund (2011) noticed 
that when is averaged at high densities, it becomes 
independent of 0o, so they took = const, which also 
provided a good fit to the density pdfs. However, they 
did not test the dependence of on M s . 

It may be more physical to retain the dependence of 
on M s even at high densities, because the postshock 
magnetic pressure tends to balance the ram pressure of 
the turbulence causing the shocks. Assuming oc A^" 1 
instead of = 0.39, we obtain eg practically indepen- 
dent of M s , as found in our AMR simulations. This 
assumption yields a density variance in MHD given by 



c x ,mhd ~ Ml /2 /2 = (T Xj hd/ M 
the ap proximation — 0.39 in |Padoan ~fc 
poTTf gave (t x ,mhd « 0.53A4 s /2 = 0.53ct x , hd 



.1/2 



for M s > 1, while 



Nordlund 
(this is- 
sue will be addressed in a future study). 

Although this correction eliminates the dependence on 
M s , the model does not predict the subtle dependence 
of eg on M a , as it does not account for the increasing 
anisotropy and dynamical alignment of magnetic and ve- 
locity fields, with decreasing M a . Therefore, the simple 
empirical law derived above is a better choice for practi- 
cal applications. 



The models by |Krumholz fe McKeel ( |2005| and |Hen 
nebelle & Chabrierl (1201 ip do not include the effect of 



magnetic fie lds. As shown by our prev ious uniform-grid 
simulations ( Padoan fe Nordlund|2011 1, even with a very 
weak mean magnetic field, corresponding to M a ~ 30, 
eg is approximately three times smaller than in the non- 
magnetized case. The AMR simulations presented here 
show that, with more realistic, larger magnetic field val- 
ues, eg decreases even further. Analytical models of star 
formation neglecting magnetic fields should not yield low 



values of eg consistent with observations of star-forming 
regions. If they do, they are inconsistent with the nu- 
merical experiments. 

5. CONCLUSIONS 

The new AMR simulations of star formation presented 
in this work show that i) eg decreases exponentially with 
increasing tg/t^ yn , but is insensitive to changes in M s (in 
the range 10 < M s < 20), for constant values of tg/td yn 
and Ma,, ii) Decreasing values of M a (increasing mag- 
netic field strength) reduce eg, but only to a point, be- 
yond which eg increases with a further decrease of M a . 
iii) For values of M a characteristic of star-forming re- 
gions, eg varies with M a by less than a factor of two. We 
propose a simple law of star formation depending only 
on iff/idyn, based on the empirical fit to the minimum 
eg: eg « e wind exp(-1.6i ff /idyn)- 

This law shows that MHD turbulence is very effective 
at slowing down the star-formation process and can ex- 
plain the low average SFR in molecular clouds. Because 
it only depends on the mean gas density and rms velocity 
of a star-forming region, the star-formation law we pro- 
pose is straightforward to implement in simulations and 
analytical models of galaxy formation and evolution. Fu- 
ture work should test its effect in simulations where the 
feedback of SNs is accounted for, but the process of star 
formation is not spatially resolved and needs to be mod- 
eled. 
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